In [1]:
# profile = 'phobos' # remote workstation
# profile = 'pantheon' # remote cluster
# profile = 'zeus' # remote workstation
profile = 'mpi' # local machine
In [2]:
import numpy as np
from zephyr.Dispatcher import SeisFDFDDispatcher
from zephyr.Parallel import CommonReducer
from IPython.parallel import Reference
In [3]:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
%matplotlib inline
import mpld3
mpld3.enable_notebook()
In [4]:
lclip = 2000
hclip = 3000
clipscale = 0.1
sms = 0.5
rms = 0.5
def plotField(u):
clip = clipscale*abs(u).max()
plt.imshow(u.real, cmap=cm.bwr, vmin=-clip, vmax=clip)
def plotModel(v):
plt.imshow(v.real, cmap=cm.jet, vmin=lclip, vmax=hclip)
def plotGeometry(geom):
srcpos = geom['src'][:,::2]
recpos = geom['rec'][:,::2]
axistemp = plt.axis()
plt.plot(srcpos[:,0], srcpos[:,1], 'kx', markersize=sms)
plt.plot(recpos[:,0], recpos[:,1], 'kv', markersize=rms)
plt.axis(axistemp)
In [5]:
cellSize = 1 # m
nx = 164 # count
nz = 264 # count
freqs = [1e2, 2e2, 3e2, 4e2] # Hz
freeSurf = [False, False, False, False] # t r b l
nPML = 32 # number of PML points
nky = 1 # number of y-directional plane-wave components
In [6]:
velocity = 2500 # m/s
vanom = 100 # m/s
density = 2700 # units of density
Q = 500 # can be inf
In [7]:
srcs = np.array([np.ones(101)*32, np.zeros(101), np.linspace(32, 232, 101)]).T
recs = np.array([np.ones(101)*132, np.zeros(101), np.linspace(32, 232, 101)]).T
nsrc = len(srcs)
nrec = len(recs)
recmode = 'fixed'
geom = {
'src': srcs,
'rec': recs,
'mode': 'fixed',
}
In [8]:
cache = False # whether to cache computed wavefields for a given source
cacheDir = '.'
parFac = 1
chunksPerWorker = 0.5 # NB: parFac * chunksPerWorker = number of source array subsets
ensembleClear = False
In [9]:
dims = (nx,nz) # tuple
rho = np.fliplr(np.ones(dims) * density)
nfreq = len(freqs) # number of frequencies
nsp = nfreq * nky # total number of 2D subproblems
cPert = np.zeros(dims)
cPert[(nx/2)-20:(nx/2)+20,(nz/2)-20:(nz/2)+20] = vanom
c = np.fliplr(np.ones(dims) * velocity)
cFlat = c.copy()
c += np.fliplr(cPert)
cTrue = c
In [10]:
# Base configuration for all subproblems
systemConfig = {
'dx': cellSize, # m
'dz': cellSize, # m
'c': c.T, # m/s
'rho': rho.T, # density
'Q': Q, # can be inf
'nx': nx, # count
'nz': nz, # count
'freeSurf': freeSurf, # t r b l
'nPML': nPML,
'geom': geom,
'cache': cache,
'cacheDir': cacheDir,
'freqs': freqs,
'nky': nky,
'parFac': parFac,
'chunksPerWorker': chunksPerWorker,
'profile': profile,
'ensembleClear': ensembleClear,
# 'MPI': False,
# 'Solver': Reference('SimPEG.SolverWrapD(scipy.sparse.linalg.splu)'),#Solver,
}
In [11]:
%%time
sp = SeisFDFDDispatcher(systemConfig)
sp.remote.dview.activate()
sp.remote.e0.activate('e0')
survey, problem = sp.spawnInterfaces()
srcs = survey.genSrc()
sp.srcs = srcs
In [12]:
%%time
sp.forward()
In [13]:
sp.forwardGraph
Out[13]:
In [14]:
%%time
dObs = sp.dPred
sp.dObs = dObs
In [15]:
%%time
sp.rebuildSystem(cFlat.T)
sp.forward()
dPred = sp.dPred
In [16]:
dResid = sp.residual
misfit = sp.misfit
In [17]:
%%time
sp.backprop()
g = sp.g
In [18]:
%pylab inline
In [19]:
ind = dResid[0].real
clip = abs(ind).max()
imshow(ind, vmin=-clip, vmax=clip, cmap=cm.bwr)
colorbar()
Out[19]:
In [21]:
ind = (g[0]+g[1]+g[2]+g[3]).real
clip = abs(ind).max()
imshow(ind, vmin=-clip, vmax=clip, cmap=cm.bwr)
colorbar()
Out[21]:
In [ ]: